#include "main.h"
#include "vec3.h"
#include "atom.h"
#include "group.h"
#include "particle.h"
#include "lattice.h"
#include "parameter.h"
//#include "gldisplay.h"


// define the global parameters
char LATT[] = "fort.15";                // file of the initial lattice 
bool PDX    =  true;                    // X Perodic b.c. condition
bool PDY    =  true;                    // Y Perodic b.c. condition
bool PDZ    =  true;                    // Z Perodic b.c. condition
int  KX     =  40;
int  KY     =  20;
int  KZ     =  40;
real M0     =  1.;                      // atom's default Mass
real XSIZE, YSIZE, ZSIZE;               // box size by lattice default
real DXI, DYI, DZI, XI2, YI2, ZI2;

// define the global groups' system
int              NUMA;                  // Number oa All particles
int              NS;                    // Number of particles' Species

vector<int>      LL;                    // Label Lower index of atoms
vector<int>      LU;                    // Label Up indew of atoms
vector<vec3>      R;                     // position vectors for all atoms
vector<vec3*>     R0;                    // all above's pointers
vector<group>    G;                     // groups for different species of atoms
vector<particle> PT;                    // Particles version of groups
sys_p            PS;                    // all MD System's Parameters

// kb-loop steps controllers 
int kb_start =  1;
int kb_end   =  10;
int kb_now   =  kb_start;

void verlet()
{}

void eval()
{}

void initialization()
{
    int i, j, ns, n=0;                  // denote total number 
    PS.initialize0();                   // initialize parameters
    lattice(LATT, NS, LL, LU, R);       // set up the lattice

    DXI = real(KX) / XSIZE;
    DYI = real(KY) / YSIZE;
    DZI = real(KZ) / ZSIZE;
    XI2 = 2.0 / XSIZE;
    YI2 = 2.0 / YSIZE;
    ZI2 = 2.0 / ZSIZE;


//  test if lattice do the right thing
    ns = LL.size();
    if(ns!=NS) cout<<"error: ns!=NS"<<endl;
    cout<<"ns is "<<NS<<endl;
    for (cout<<"LL is ",i=0; i<ns; i++) {
        cout<<LL[i]<<" ";
    }
    for (cout<<endl<<"LU is ",i=0; i<ns; i++) {
        cout<<LU[i]<<" ";
    }
    NUMA = R.size();
    cout<<"\n the size of R is "<<R.size()<<"\n";
    cout.precision(2);
    cout.setf(ios::fixed);
    cout<<"position of first atom: "<<R[1].x<<" "<<R[1].y<<" "<<R[1].z<<endl;
    cout<<"position of last  atom: "<<R.back().x<<" "<<R.back().y<<" "<<R.back().z<<endl;

    G.clear();
    G.resize(NS);
 
    for(i = 0; i < NS; i++){
        G[i] = group(i+1, LL[i], LU[i], 1);           // generate groups
        G[i].r1_random();
        //G[i].info();
    }

   R0.push_back(&VO);
   for(i = 0; i < NS; i++){
        for(j = 0; j < G[i].num; j++){
            //G[i].a[j] = atom(i+1, LL[i]+j, M0, R[n]); // set up positions
            //R0.push_back(&(G[i].a[j].r0));            // and their pointers
            R0.push_back(G[i].a[j].getR0());            // and their pointers
            n++;                                        // total number++
        }
    }


    n = 0;
    PT.clear();
    PT.resize(NS);

    for(i = NS-3; i < NS; i++){
        PT[i] = particle(i+1, LL[i], LU[i], 1);        // generate groups
        for(j = 0; j < PT[i].num; j++){
            PT[i].a[j] = atom(i+1, LL[i]+j, M0, R[n]); // set up positions
            //R0.push_back(&(PT[i].a[j].r0));          // and their pointers
            n++;                                       // total number++
        }
    }


    for(i = NS-3; i< NS; i++) PT[i].info(); 


    // to test if R0 do the right thing
    if(R.size()!=R0.size())cout<<"error: R0 size ="<< R0.size()<<endl;    
    for(i=0; i<NUMA; i++){
        if((R[i]-*R0[i]).module()>0.00001){
            cout<<"error: R[i] != R0[i] "<<i<<endl;
            break;
        }
    }
    for(i = 0; i < NS; i++)
        for(j = 0; j < G[i].num; j++)
        {
            vec3 ri = *(G[i].a[j].getR0());
            vec3 rj = *R0[LL[i]+j];
            if((ri-rj).module()>0.0001)
            {
                cout<<"error: R[i] != R0[i]  at i="<<LL[i]+j<<endl;
                exit(3);
            }
        }

    G[1].a[24].info();

    /*
    cout<<"R.size()  = "<<R.size()<<endl;
    cout<<"R0.size() = "<<R0.size()<<"\n";
    cout.precision(2);
    cout.setf(ios::fixed);
    cout<<"position of first atom: "<<R[1].x<<" "<<R[1].y<<" "<<R[1].z<<endl;
    cout<<"position of last  atom: "<<R.back().x<<" "<<R.back().y<<" "<<R.back().z<<endl;
    cout<<"position of first atom: "<<(*(R0[1])).x<<" "<<(*(R0[1])).y<<" "<<(*(R0[1])).z<<endl;
    cout<<"position of last  atom: "<<(*R0.back()).x<<" "<<(*R0.back()).y<<" "<<(*R0.back()).z<<endl;
    // generate pgroup
    // pgroup pg()
    */
}

void kb_loop()
{

    int i;
    
    for(kb_now = kb_start; kb_now <= kb_end; kb_now++)
    {
        for(i=0; i<NS; i++) G[i].predict();
        eval();
        for(i=0; i<NS; i++) 
        {
            G[i].correct();
            G[i].periodicbc();
        }
        G[1].a[24].r0.display();
        
    
    }

}


int main(int argc, char **argv)
{

    initialization();
    //kb_loop();

    return 1;
}

